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ABSTRACT 


The optimal minimum time control (1.e. bang-bang controller) is applied to 
the fast reaction missile defense problem. From Pontryagin, the optimal control 
was determined to be a function of the adjoint in the minimization of the 
Hamiltonian. The control may also be posed either as a function of time or as a 
function of the states. The state space can be partitioned into regions, surfaces 
and curves where the optimal control action is either its maximum plus or minus 
N. 

In missile simulation problems, the method of adjoints 1s often used in 
parametric studies of errors and miss distance. This technique is developed 
both graphically and mathematically, and is used here to help one visualize the 
solution trajectory and families of optimal trajectories for all possible initial 
conditions. 
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I. INTRODUCTION 


In this era of technological advancement, more and more effort is being 
devoted to automation to increase speed and efficiency. Everything is 
becoming faster, smaller, more efficient and more effective, from fuzzy logic 
controlled appliances to smart weapon systems. As our engineering and design 
of physical systems improves we are able to generate faster and more accurate 
mechanical devices, and the control systems for such devices must improve as 
well. In pushing to develop the fastest, most efficient controller for whatever 
our application, we introduce the bang-bang controller. 

Many of our current defensive missile systems were designed to be used 
against threats that are no longer of greatest Importance. A point defense 
missile may now be required to bring down a target that has a significant speed 
advantage, and be able to do it in such a way as to control the fallout after the 
impact. 

Even if we redesign our systems, current economic conditions make it 
unlikely that we could build up an arsenal of new weapon systems. However, 
we could upgrade our currently existing arsenal by replacing the original control 
logic with a newer, more capable controller. 

In this report we will develop second and third order minimum time optimal 
controllers and apply them to a fast reaction missile defense problem. While we 
are using the missile defense problem as our example, it should be noted that 
the controllers are not limited to this example and may be applied as generic 


controllers for a wide variety of systems. 


II. SECOND ORDER CONTROLLER 


A. MINIMUM TIME OPTIMAL CONTROL 

Designing a control system is frequently a hit-and-miss process using a 
variety of design techniques to iteratively create a system that meets specific 
criteria. The performance of the system may be defined in the time and 
frequency domain in terms of overshoot, settling time, etc., or it may be defined 
by some externally measured criteria. For example, a satellite control system 
may be designed to minimize fuel consumption. 

"The objective of optimal control theory is to determine the control signals 
that will cause a process to satisfy the physical constraints and at the same 
time minimize (or maximize) some performance criterion." [1] 

1. Problem Definition 

The system and optimization problem is defined as 
Хх = Ах + Ви (2.1) 


with x(O) =c, and x(t¢) = 0, while minimizing 


J= fat (2:2) 
0 


Let us consider the simple example with 


X - ТЫН o 


From Pontryagin [1], we find we can minimize J by minimizing the Hamiltonian 


IT— T Opp (2.4) 


This is minimum when u is operating at its maximum possible value and with 


opposite the sign of the adjoint ps. Thus we have 


u--N-sign(p5) (2.5) 
where 
p=- ol (2.6) 
This has a solution 
ри (255) 
Pp =—Cyt+Cp. (2.8) 


A typical solution would appear as seen in Fig. 2.1. 





Figure 2.1 Solution to Minimized Hamiltonian 


Note that from the adjoint solution, the control may change sign only 
once. We would like to solve this problem for all possible initial conditions and 
only one terminal condition (x(ts) = 0). Hence it makes sense to look at this 
problem in negative time (adjoint), starting from the end point of motion. From 


the uniqueness theorem in ordinary differential equations, only two trajectories 


can emanate from the origin; one for u = -N and one for u = +N, as shown in 


Figure 2.2. In this second order example, our solution is constrained to the x, 


х2 plane. These negative time trajectories divide the state space (here a plane) 


into two parts. Solution trajectories emanating from these curves constitute all 


possible solutions for all possible initial conditions. 





Figure 2.2 Zero Trajectory Curves 


2. Minimum Time Trajectories, Second Order 


ЫК 


This can be represented in flow diagram form as 


Consider the system 


(2.9) 


The control, u, may be defined for a bang-bang control system as +N, where N 
IS some constant. Given any starting values for the states, there are only two 
possible paths for the states to take; one corresponding to u = +N, and опе 
corresponding to u = -N. If we desire to drive all states to zero, and we know 
that we can only apply u = +N, Figure 2.3 shows the paths followed by the 


states given various starting values. 





Figure 2.3 Minimum Time Trajectory Solution Curves 


3. Second Order Solution Trajectories 
This system has a solution 
x(t) = $(t,0)x(0) -- A(t,O)u(O) (2.10) 
where 
ф=е“, (2.11) 


and 


A = [ева (2.12) 


Expanding e^', we get 


2! 


pa Го а 

= + t+— t +... 
0 1 О) 2!|0 00 0 
ОЕ l t 

= E m: t +... = | (22059 
О 1 ОТОО Ore 


Since A? is the zero matrix, all higher orders of A will also be the zero matrix 


c2 s(a A+.. 


and we may drop them. Evaluating A (for u fixed) we find 


1 t—7t | 0 
۸ CCT dt 
0 A 


t-1 ue [42 
| je- P s?) (2.14) 
0| ] T Р t 


The solution for this system is 


: j Сі 
x(t)= n SI qn (2.15) 
QO ] t 
or in scalar equations 
x, (t) 2 x, (0) x4 (0)t- 3u(0)t* (2.16) 
x4 (t)  x4(0)-- u(O)t. (2.17) 


These equations describe the states as a function of time given any initial 


conditions and the fixed control effort, u. 


B. ANALYTIC SOLUTION FOR SECOND ORDER SWITCHING 
TIME 


We may now treat this system as a boundary value problem and 


analytically solve for switching times of the control effort. 


1. Solving Boundary Value Problems 

Since the control is piecewise constant, (+N), we can Separate the 
problem into two pieces and match boundary values at the point where u 
changes sign. In other words, if the system moves from point x(to) through 
point x(t,) to point x(t), we can solve the problem in two parts. Each of our 
boundary value problems can be stated in such a way as to supply simplifying 
boundary conditions, i.e., setting initial or final values to zero. Optimal control 
theory tells us that for a second order system there will be at most one 
switching time, the change from u = +N to u = -N, or visa versa. Let us 
consider first the solution for our system with some initial condition [x4(0), 


x2(0)], as shown in Figure 2.4. 


x (9), x0) 





Figure 2.4 Minimum Time Trajectory From a Fixed Initial 
Condition 


The time from t=0 to t, is the length of time before the control effort 
changes sign. We cannot solve directly for t, because we do not have enough 


information on the boundary values. Since we are solving this problem for 


arbitrary initial conditions, x1(0) and x4(0), x4(t;) and x2(t;) are unknown. We 
do know, however, that the final states, x4(tr) and x(ts), will be zero, and from 
this we can determine the zero-trajectory curve in negative time from the 
origin. This curve, shown in Figure 2.5., is the only curve that will go through 
the origin with u = +N. Therefore the switching time will occur when the path 
of a u = -N curve intersects this curve. There are an infinite number of u = -N 
curves intersecting this curve, however, only one curve will go through any 
particular set of initial conditions. 

To simplify the problem, we will start with initial conditions on the x, 
axis. Define x,(Q) = c), x2(O) = O, x,(ts) = xo(ts) = O, and u=-N. This situation 
is described in Figure 2.6 where t; may be defined as the time from tg to ty. 


Equation (2.16) may now be written 


x (t) =x,(0)+x,(0)t+4ut* = x,(0)-4Nt? (2.18) 
Since the curves for u = +N and u = -N are symmetric, it can be noted that 
xi(t) 7 3x,(0) (2.19) 
and therefore 
Lx (0)=x,(0)-+Nt/° (2.20) 


and finally 


а үз”. CA 


This is valid for any initial conditions such that xi(0) > 0, x2(0) = 0 and u = -N. 





Figure 2.5 Simplified Boundary Value Problem #1 





Figure 2.6 Simplified Boundary Value Problem #2 


The next simplification involves getting from the positive x2 axis down 
to the positive x; axis as shown in Figure 2.7. Given x4(0) = 0, x2(0) = positive 
real, x (tr) = positive real, x(t) = 0, and u = -N. 


Equation (2.17) may be written 


х,(1) = Ое (2.22) 
Аїї= 6, х(6) = 0 giving 
O X, (0)— NG (2.23) 
and 
х,(0) 
= = 2.24 
2 = (2.24) 


The final stage is to translate the initial condition off of the x2 axis to 
some unknown initial condition. The time required for the states to get from 
some initial condition x4(0) to x»(t) is not based on the x, state, and is therefore 
always equal to t», as seen in Figure 2.8. 


We solve (2.16) for х (1) 


xi(t5) =x,(0)+x,(0)t, +} Nt; 


0 0)Y 
x0) 022) -pn 20) 


| | 
= x1 (0) + х2 (0). (2.25) 


10 





Figure 2.7 Simplified Boundary Value Problem #3 





Figure 2.8 Simplified Boundary Value Problem #4 
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Finally we combine these solutions to determine the switching time, ts, 
for any initial conditions such that x4(0) » 0 and x4(0) » O with u 2 -N. From 
Figure 2.3 we see that x,(O) for the tı solution is the same point aS x;(t2) from 


the t? solution so that 
X CO) ECT) 


VN N 





l a 
коо 0) 


М 
= i (2.26) 
N AFN 
and 
t =t +t, = 0), 1 ١ E (2.27) 
N 2\ N N 


These equations were derived for specific initial conditions and are 
only valid where the initial conditions lie above the zero trajectory curves of 
Figure 2.2. In order to expand the capabilities of the control system for all initial 
conditions we re-derive the switching time equations for initial conditions 


x1(0) «€ 0, x5 (0) « O, and u(0)= +N for the solution 


в 0) 16409) o (2.28) 
N 2X N N 


This equation is valid only when the initial values are below the zero trajectory 
curves of Figure 2.2. For this discussion we will limit ourselves to initial 


conditions above the zero trajectory curves and use the solution (2.27). 


12 


2. Simulation of Analytic Switching Time Controller 
Using a computer simulation, we test the accuracy of the solution by 
choosing initial conditions and observing the response to our control effort. We 
simulated the example system of (2.9) with x4(0) » 0, x2(0) » 0, and N = -1. 
The control effort, ип, changes sign to -N at the calculated switching time, ty. 
The results show that the states pass through the origin of the state space and 
the error at the origin 1s associated with the discretization of the simulation, as 
shown in Figures 2.9 and 2.10. The switching time shown in Figure 2.10 is the 
calculated value. The control effort actually switches at the first sampling time 
following the caluculated switching time, t,. If we increase the sample rate for 
the simulation we reduce the terminal miss error. Upon reaching the origin 

some additional control logic must be devised to maintain position. 
While the switching time solution is a minimum time solution for driving 
the states to the origin, or any desired values, it must be shut off at tr. The 
switching time solution cannot adapt to a time varying situation as the switching 


time is defined solely on the initial conditions. 


Phase plot - Switching Time 
x, x, (0), x2(0) 





Figure 2.9 Simulation of Analytic Switching Time Solution 


Control Effort - Switching Time 


Шш МЕ 


Типе ($ес) 





Figure 2.10 Control Effort for Analytic Switching Time 
Simulation 
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C. SOLUTION TO SECOND ORDER SWITCHING CURVES 

An option for improving the capabilities of the second order controller 1s to 
remove the time dependancy of the control effort. The control problem may be 
posed as a function of the states Instead of a function of time for some simple 
systems. By defining the control effort as only a function of the states the 
control effort is updated continuously and may change immediately in response 
to a change in the system parameters. 

There are only two possible control efforts, -N and +N. Therefore, any 
position in space will have either +N or -N control effort to drive it along its 
unique path to the origin, and we can solve this system for the second order 


Switching curves that divide the state space into two parts, see Figure 2.2. 


1. Second Order Switching Laws 


Given the system 


Xi 0 1| x 0 
E +| и (2.29) 
X5 0 0| x> 1 
we can integrate the state equations. Setting u = +N 
x(t) =] x_(t)dt = | udt 2 £Nt «c, (2.30) 
x (t) = ОШ = | ха E | £Nt * c, di = КМ To De (2.31) 


Choosing the boundary values so that xj(ts) = x2(tp) = O we may rewrite the 


equations as 


X (t,)= + ЕМИ + (0) x, (0) (232) 
X5 (tp) 2 €ENt, t x5 (0) (2:931) 

Or 
О= +1 №? +х,(0); +х,(0) (2.34) 


O=+Nt + x» (0). (2.35) 


Solving for tr 























_ Хо (0) 
ЕЕ T (2.36) 
Then for any value of t 
L (2.37) 
N 
Substituting this into equation (2.34) 
х (1) x(0 Y 
0= 00+ 0019 2 ЈЕ : | (2.38) 
N N 
2 2 
— X3 t) x5 (t) 
xi (t) N 2 N ( ) 
x3 (t) 
O=x, (t) ++ — —. (2.40) 
N 
We make the substitution 
+x2(t)=-x,(t)|x>(t)| (2.41) 
so that 
0 = x, (t)- şk x; CO|xs (t)]. (2.42) 


Equation (2.42) describes the two parabolic trajectories referred to as 
the zero traJectory curves because any states on these curves will, with the 
appropriately signed control effort, be driven to the origin. These zero 
trajectory curves divide the state space into two regions of opposite control 
effort. If the states are above the curves then u = -N, but if they are below the 


curves then u=+N. Therefore our control effort, u, may be defined 


u =—N-sign(x, (t)— 5k xo CO|xo (0|). (2.43) 
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With this switching law we drive the states from any initial conditions to the 
origin with at most one change of the control effort. 

For example, let us define the initial conditions of the states to be 
above the zero trajectory curves so that x4(0) « O and x,(0) > 0, as shown in 
Figure 2.11. From (2.43) the control effort is u = -N, which drives the states 
along the parabolic trajectory shown in Figure 2.11. 

When the states intersect the zero trajectory curve, the control effort 
changes, according to (2.43) to u = +N, and follows the zero trajectory curve 
into the origin. 

2. Limit Cycles 

The control effort, u = +N will not only drive the states to the origin, it 
will, in fact, become confused there. It is a point of indecision and chatter or 
limit cycle motion ensues as in Figure 2.12. The magnitude of the limit cycle 
depends on the sampling rate or time delay for the control effort. If the sample 
rate is low, the states will penetrate far into the opposite control region before 
the control effort can change, and the limit cycle will swing widely around the 
origin. If the sample rate is high then the states will not travel far away from 
the origin before the control effort corrects the direction of travel. 

Since a heavy chatter mode is not desirable for many real systems, 
control logic for reducing or removing the chatter mode may be developed, 


however, it will not be included in this paper. 
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x1(0), x; (0) 


-N 





Figure 2.11 Solution Trajectory From Fixed Initial Condition 





Figure 2.12 Limit Cycle on Second Order Solution Trajectory 


3. Simulation of the Switching Law Controller 
To test the switching law (2.43), we simulate the system of (2.29) 
using a maximum control effort of N = 1 and initial conditions x,(0)>0, 
x,(0)>0. The output of the simulation, shown in Figure 2.13, demonstrates the 
control effort driving the states first to the zero trajectory curve, then along the 
zero trajectory curve to the origin. Once at the origin, the control effort goes 


into limit cycle as shown in Figure 2.14. 


КО 


Phase plot - Switching Law 


. х (0), х (0) 


Possible Chatter 
or Limit Cycle 





Fieure 2.13 Simulation Using Second Order Switching Law 


Time (sec) 





Figure 2.14 Control Effort for Second Order Switching Law 
Solution 


III. APPLICATION OF A SECOND ORDER CONTROLLER 


A. MISSILE/TARGET INTERCEPT MODEL 

An application for Minimum Time Optimal Control is an anti-missile defense 
system, where the incoming target has a speed advantage over the defensive 
missile. In this case the missile control system must operate in saturation mode. 
The bang-bang controller, where control effort is either maximum-positive or 
maximum-negative, has a faster response capability than the standard 
Proportional Navigation guidance system. A simple model of a missile to target 
engagement can be described by Figure 3.1 where o is the line of sight (LOS) 
angle from missile to target, y, 1s the angle of the missile velocity, yı is the angle 


of the target velocity, and all angles are relative to an inertial reference. 


Missile 





Figure 3.1  Missile/Target Geometry 


It is assumed that the target is on the final leg of it's flight and is now on a 
straight, non-maneuvering traJectory. The control is to drive the line of sight 
гаје (6) and it's derivative (G) to zero in minimum time. 

For our models we will be using only two dimensional scenarios. The 
system dynamics use second order models for each dimension. The target 
dynamics are non-maneuvering so that the acceleration, a, is zero. The missile 


acceleration, aj, is assumed perpendicular to the velocity vector yg. The 


system dynamics are shown in the signal flow graph in Figure 3.2. 


o 
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Figure 3.2 System Dynamics of the Intercept Model 
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The LOS angle, o, and it's derivatives, 6 and 6, may be derived analytically 
from the available states. We define the geometry of the system as in Figure 


B 


Missile 





Figure 3.3 Geometry for Angle Definitions 


with 
R, = Range of the Target. 
Ка = Range of the Missile. 
У, = Velocity of the Target. 
Ма = Velocity of the Missile. 
a, = Acceleration of the Target. 
да = Acceleration of the Missile. 
Ve = Closing Velocity. 
so that 


R 
O = ал | к" | (3.1) 


tmx 


R У - К, У 
СЕ unx "tmy | tmy ‘tmx (3.2) 
R 
and 
.. Rmx аш У c nn ua Па Ven Rae oe 
О = Шо فك‎ a ЕЕЕ 3 (3.3) 
К К 
where 
М, =-Е (3.4) 


In an actual application Kalman or Luenberger Observers may be used to 


obtain estimates of 6 and 6. 


B. APPLICATION OF THE SWITCHING LAW 
The switching law from (2.43) 1s applied to our scenario where б апа б 
form the state space of this simulation, and is implemented as 
: TWO uo 
ie зна e ل‎ | (3.5) 
The sign convention of (2.43) was changed to conform to the geometry of our 
Scenario. 


We set the simulation with the initial conditions 


Ко 00) = О Rmy(O) = 0 ft 

У (0) = 2000 ft/sec У пу()) = 0 ft/sec 
a, @ (0) = 6С ату(0) = 0 ft/sec2 
R,,(0) = 10000 ft R,y(O) = 1000 ft 
V, (O) = 20001 see V, ,(0) = 0 ft/sec 
3 (O) = Oil) sec a, (0) = 0 ft/sec? 
i12 29 Sec dt = 0.01 sec 
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Running the simulation we find that the controller does drive 6 and 6 to 
zero, and maintains them about the origin until intercept is reached (Figure 3.4). 
Once the states reach the origin, the system controller goes into chatter mode, 
see Figure 3.5, and the trajectory chatters back and forth about the desired 


intercept path. 
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Figure 3.4 





Figure 3.5 
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Order Controller 
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Control Effort for Second Order Missile/Target 
Simulation 


26 


C. COMPARISON WITH PROPORTIONAL NAVIGATION 
CONTROLLER 


The system dynamics for the proportional navigation controller simulation is 
basically the same as that in the previous section, except that the control effort , 
u, is proportional to 6, which is obtained with an estimator, shown in Figure 3.6, 
where B =ó. We limit the acceleration so that u is bounded by +N. With only 
the change in the controller we ran the simulation for the same initial conditions 
with the results shown in Figures 3.7 and 3.8. 

This Proportional Navigation system is effective at intercepting the target 
only when the missile/target geometry and kinematics are sufficient that the 
missile has time to maneuver. The delay in saturation of the control effort 
caused the missile to turn too slowly so that the Proportional Navigation 
Controlled system was unable to maneuver quickly enough to intercept the 


faster, although non-maneuvering, target. 
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Figure 3.6 System Dynamics of Proportional Navigation 
Missile/Target Model 
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Control: Proportional Navigation 


Intercept Time = 2.02 sec 
Miss distance = 370.3 ft 
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Figure 3.7 Proportional Navigation Missile/Target Simulation 
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Figure 3.8 Control Effort for Proportional Navigation 
Missile/Target Simulation 


IV. THIRD ORDER CONTROLLER 


A. FORWARD TIME SYSTEM 

The second order controller is effective at intercepting a target that has 
speed advantage. Previously we have controlled the system by driving © and 
O to zero, and thereby maintaining a constant LOS angle, o, until impact. But 
what about controlling the LOS angle itself? There are situations in which we 
would like to be able to attack a target from a particular angle, or perhaps have 
multiple missiles, each with it's own pre-defined attack angle. 

Consider a point defence system in which the missile, launched from some 
point away from the target's final trajectory, would try to position itself in 
minimum time onto a "head-on" collision course with the target, as shown in 
Figure 4.1. In such a situation we would use a third order minimum time 
controller to drive б, б, апа G to zero. 

1. System Definition 


Our third order example is 


0 1 0 0 
хе=| О ONUS (4.1) 
0 0 0 1 


Minimizing the Hamiltonian and solving the system we find 
Н =1+рх, + рх, +р;и (4.2) 
апа 


у = —N-sign(p;) (4.3) 
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Attacking Missile 
Trajectory | = 
Attacking Missile 


Target Ship Intercept Point 


Defensive Missile 


очи Minimum Time solution for 


a zero LOS intercept angle. 


Escort Ship 


Figure 4.1 Application of Third Order Minimum Time Control 
where 
€ 0 0 0 
р=-—=10 -1 0 |р 
Ox 
0 0 -1 
so that 
Pie 


рә = —Cit + C5 


p, 2 ict -ct6,. 
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ler 


(4.4) 


(4.5) 
(4.6) 
(4.7) 


From the adjoint solution the control may switch sign no more than 
twice. Again tracing this problem in negative time we may follow the zero 
trajectory curves out from the origin with control efforts of £N. These curves 
are now in three dimensional space residing in their own plane which intersects 
the origin. Intersecting these zero trajectory curves are an infinite number of 
curves making a surface, and leading off from this surface the infinite number of 
trajectories lead to the initial conditions. Therefore in forward time we may 
start with an initial condition such that u = +N will drive the system to intersect 
with the surface at tsy as shown in Figure 4.2. Switching the control effort to 
u = -N will drive the system along the surface to intersect with the zero 
trajectory curve at t.w2 where the control effort switches again. Finally u = +N 


will drive the system to the origin. 


Хү 


Zero Trajectory 
Curve 


Curved 
Surface 





Figure 4.2 Three Dimensional Switching Curves 


2. Third Order Switching Curves 
We have determined from the second order system that solving for the 
control switching times is not as widely applicable as defining the control as a 
function of the states; so we now move on to defining the third order switching 
curves. 


Starting with the state equations (4.1) we discretize and expand the 


equations so that 





NO 253 
ES EE 
! 3! 
и ао 0 1 0 0 0 1 
Е ООО РОО Е ООО 
OO 0.0 0 ОО 
1 t dic 
SrO i t (4.8) 
0 0 1 
A? is the zero matrix so it and all higher terms of A are dropped. 
| | t-t 1(ї-т)°|[0 
Az | ем“ ват = | О 1 t-1 | О [ат 
0 0 
() ] ] 
t 
E e 
=| t— due t-il = 117 (4.9) 
] Т t 
0 


Combining these we obtain 
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2 3 


l t 5t +i 
x()=|0 1 t |x(0)+| 1 |и(0) (4.10) 
0 0 1 t 
or In scalar equations 
xi (t) 2 xy (0) tx (0) - 3 6x4 (0) - 1 6u(0) (4.11) 
x) (t) 2 x4 (0) tx4(0)-- 11^u(0) (4.12) 
X4 (t) 7 x4(0)- tu(0). (4.13) 


B. NEGATIVE TIME SYSTEM 

The third order system, being piecewise continuous, is easily broken into 
several simple boundary value problems. In order to solve the systems of 
curves we must determine some boundary values for the intersections of these 
families of curves. We start at the origin and run the system in negative time to 
determine our other boundary conditions. 


In negative time we have 


0 -1 0 0 
x= O O -1x-|O yu. (4.14) 
0 0 0 -] | 
Discretizing we find 
l — 141° 
þ=e^™=|0 1 -t (4.15) 
0 0 1 


and 
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1 
t 
A= | eX Badr = 12 (4.16) 
0 
= 
so that 
2 3 
1 -t 4t -4t 
хи)=|0 1 = |х(0)+ 21 |00), (4.17) 
0 0 ] -t 
or in scalar equations 
xı (t)= x (0)- tx, (0)+4t°x;(0)-4tu(0) (4.18) 
x) (t) 2 x,(0)—tx,(0)+417u(0) (4.19) 
X4 (t) 2 x4(0) - tu(0). (4.20) 


1. Solving for Negative Time Boundary Points 
To develop a complete solution we solve the equations for several 
different points along the zero trajectory curve. Setting u = -N and traveling out 


along the zero trajectory curve for 1 second we find 


x (1) 2 - 1uc 2-1(-N)2 iN (4.21) 
х,(1) =1ш? =1(-№)=-1М (4.22) 
хз(1) = -ш= -(-№ = №. (4.23) 


Similarly we run the system in negative time for 2 and 3 seconds to 


obtain other boundary points as shown in Figure 4.3, and listed in TABLE 1. 


EN 


Zero Trajectory 
Curve 





Figure 4.3 Boundary Values in Negative Time Solution 


NEGATIVE TIME BOUNDARY CONDITIONS, u = -N 





TABLE 1. 


The control effort may also be u = +N, traveling away from the origin 


on the other zero trajectory curve giving the boundary conditions in TABLE 2. 
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NEGATIVE TIME BOUNDARY CONDITIONS, u = +N 





TABLE 2. 


We may now solve the equations for the family of curves that intersect 


the zero trajectory curves. Specifically, we solve for the equation of the one 


curve that travels from some point x4(0), x2(0), x3(0), to the point x4(t) — iN, 
x(t) = +N, X3(t) = N. Since the control effort on the zero trajectory curve for 
this intersection point is u = -N, the control effort for the curve we are solving 


for must be u = +N, and the forward time equations may be written as 


xı (t)=}N =x (0)+tx,(0)+4 tx, (0)++ tN (4.24) 
x. (t)=-+N=x,(0)+tx,(0)++N (4.25) 
хз(1) = М = х. (0) +103 (4.26) 


Solving (4.26) for t we find 


QUE 


t (4.27) 


Substituting (4.27) into (4.24) and (4.25), and after a bit of algebra we obtain 


2 
-N = x2(0)- ŽA (4.28) 
2 3 
ТО со л, (429) 


N 2N  3N° 
Using the other boundary values from TABLE 1 and TABLE 2 as 


well, we may generate a family of equations representing both sides of the zero 


ЭЙ 


traJectory curves, and different points of Intersection along the curves (see 
TABLE 3). The control effort, u, in TABLE 3. 15 the control effort for the zero 
trajectory curve that is Intercepted. The time, t, is the time out from the origin 
to the intercept point for the negative time system. 


To combine all the equations from TABLE 3 into one solution we 





define 
В ха [хз 
w = slen| x, + 4.30 
£ | 2 2N | ( ) 
апа 
а 
f = х, + №. 4.3] 
2 2N s», 


so that our family of curves is defined as 

0= ×, )0(+ Qu s Qr QD, t TE (4.32) 

N N 

Equation (4.30) determines which signs are to be applied based on which 
direction the zero trajectory curve is on; Equation (4.31) adjusts the magnitudes 
of the equation depending on the distance of the intersection of the zero 
trajectory curve from the origin. Each curve intersecting the zero trajectory 
curve will have a different value of the function f, and f will remain constant all 


along that curve. 
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FAMILY OF SOLUTION CURVES 
2 
-N =x2(0)- ха (0) 
2N 


2 3 
O=x,(0)+x,(0) x5(0)x4(0) X3 (0) | X3 (0) 


Хз ИО) * (0) 
aN 3 

O=x,(0)+2x,(0) x5(0)x4(0) X3 OF UA 
x (O) 

2N 
0 2 x1(0) - 3x4 (0) 
O 

2N 
02 x,(0)-x,(0)— 


2. "> 


-4N - x2(0)- 


—ON 2 x2(0)- 
x2(0)x (0) 3x47 (0) | x; (0) 


N = x2(0)+ 


x2(0)x (0) x > (0) x (0) 
N 2N 3N? 


4N = x2(0)+ `—— 


x2 (0)x3(0) | x (D) E x4 (0) 
N N 3N? 


О=х,(0)+2х„(0)+ 


х (0) 
2 N 


О=х,(0)+3х„(0)+ 


QN = x2(0)+ ——— 
x4(0)x4(0) E О) И (5) 

N 2N 3N? 
TABLE 3. 





C. THIRD ORDER SWITCHING LAW 

Using (4.30), (4.31) and (4.32) we can break up space into areas of 
opposite control effort. The control effort is defined by which of these volumes 
contain the states. Therefore, the third order switching law is defined as 


0 0 0 
u=-N: sies 582. ЈА „Меле, E) (4.33) 


D. THIRD ORDER CONTROLLER SIMULATION 

A simulation of this third order model shows that the third order switching 
law, (4.33), drives the states to the origin with only 2 changes is the control 
effort. In application some means of shutting off the control effort will be 
required for the system to avoid entering chatter mode upon reaching the origin 
(see Figure 4.4). 

An examination of the functions (4.30), (4.31), and (4.32) in Figure 4.5 
shows that f is a smoothly increasing curve until the states intercept the 
switching surface. Once the states are following the surface, f remains a 
constant value. When the states reach the zero trajectory curve, f becomes 


approximately zero. 
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Figure 4.4 Simulation of a Third Order Minimum Time 
Controller From a Fixed Initial Condition 
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Figure 4.5 Control Laws For Third Order Solution 
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V. MISSILE ADJOINTS 


Given a time varying system 


x — A(t)x * B(t)u (5219 
у = С(0)х. (5.2) 
и сап be shown that the impulse response of the adjoint system 
Be A(t pene ela (5.3) 
T 
y=B(tr—t) p (5.4) 


is the response of the original system, at time t, to an impulse applied at time 
t; —t before t. [2] For a time invariant system the transfer function for the 
original system is identical to the transfer function of the adjoint system, Le., 


they are self-adjoint. 


A. CONSTRUCTION OF AN ADJOINT 

Given the mathematics of the adjoint method, we now need to define some 
rules to create and understand the adjoint system with real systems and block 
diagrams. [3] 

1. Rule 1: Convert All System Inputs to Impulses 

In constructing the adjoint it is necessary that all system inputs be 

impulsive. Since this may not be the case with the block diagram, all system 
inputs and initial conditions must be converted to impulsive inputs via block 
manipulations and extensive use of integrators.. Figure 5.1 shows that step 
inputs and initial conditions are equivalent to the output of impulse driven 


integrators. 


2. Rule 2: Replace t With tí - t in the Arguments of All Time 
Varying Coefficients 
In many linear systems it is possible to express a gain as a function of 
time. The adjoint system operates in negative time and Figure 5.2 shows the 


conversion of functions of time into the adjoint domain. 





Figure 5.1 Conversion to Impulsive Inputs 


Original System Adjoint System 


Time K(t)=at+b K(t—t;) =a(t—t,)+b 


Varying 


Gain l 
K(t)=—— К(ї—ї)= 
s b te) at +b 





Figure 5.2 Convert Functions of Time to the Adjoint Domain 
3. Кише 3: Reverse the Direction of all Signal Flow 
The direction of all signal flow must be reversed, redefining nodes as 


summing junctions and visa versa. Notice that all system outputs become 
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inputs and all system inputs become outputs. This last rule allows for the simple 
graphical creation of the adjoint system by first drawing the block, or flow, 
diagram of the original system, then redrawing it with the arrows reversed. 
Figure 5.3 shows some examples f converting nodes to summing junctions and 


visa versa. 


Original System Adjoint System 





Figure 5.3 Converting Nodes to Summing Junctions 


B. DEVELOPMENT OF A SECOND ORDER ADJOINT SYSTEM 

This adjoint formulation lends itself well to analyzing some optimization 
problems, those where tr is free and the terminal state is constrained to a point, 
curve or surface. From Pontryagin the control is a function of the homogeneous 
adjoint (see (4.2) and (4.3)). The missile adjoint form gives also a particular 
solution of the adjoint (1.e. system impulse response). 

The adjoint allows one to generate optimum solutions (i.e. switching 
surfaces) for all possible initial conditions. The forcing impulses passed through 
an integrator gives our saturation type of optimal control (XN). 

As an example for the application of the method of adjoints we will apply 
our adjoint rules to our second order controller. The system described by (2.9) 


has a time dependent control input where 
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NEUES E. Be 
OS +N t>t, d 


and, from our second order example, 


m. 0) 1020), NE (5.6) 
N 2l N N 


1. Apply Rule 1 to the Second Order Example 
For this example we will define N = 1. We first draw out the original 
signal flow diagram, showing all the inputs and initial conditions. In order to 
apply the impulsive inputs we expand the second order system to the third 


order state equations 


010 0 
X210 0 I|x4|O lu (5.7) 
000 ] 
у=[1 0 O]x. (5.8) 


Converting all the inputs to impulsive inputs and initial conditions we get the 


flow diagram in Figure 5.4. 
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Original System Showing All Inputs 
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Figure 5.4 Application of Rule 1 on Second Order Example 
2. Apply Rule 2 to the Second Order Example 


We replace t with t - t; in the arguments of all time varying coefficients 


to obtain the flow diagram in Figure 5.5. 
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oi Jc O UC 
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Figure 5.5 Application of Rule 2 on Second Order Example 
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3. Apply Rule 3 to the Second Order Example 
Finally we reverse the direction of the signal flow, redefining nodes as 
summing junctions and visa versa, thereby changing the system inputs to 


outputs and the system outputs to inputs as shown in Figure 5.6. 


Adjoint System 


p /uwt-t)i m A 000—0) (0) 





Figure 5.6 Application of Rule 3 on Second Order Example 
The mathematical definition of the Adjoint System is 
p=A'p+Clu (5.9) 
y=B'p. (5.10) 
Using the A, B, and C matrix from (5.7) and (5.8) we obtain 


Pi 00 Ор: l 
p |= 1 O O| p. I+ O [u(t—tk) (5.11) 
P3 0 1 0 P3 0 

у=[0 0 1]р (5.12) 


which corresponds to Figure 5.6. 


C. SIMULATION OF THE SECOND ORDER ADJOINT 
1. Forward Time Second Order Simulation 
In the forward time second order simulation we set N = 1, and the 


initial conditions x(0)=2 and x(0)=0. Using impulses through integrators we 
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apply the initial control effort at t = 0. A second impulse, opposite in sign and 
twice in magnitude, is applied at the switching time, as defined in (5.6), driving 
the states through the origin. 

The simulation length is 4 seconds with a sample step of .05 seconds. 
The results are shown in Figure 5.7. Starting at the initial conditions the system 
is driven through the origin with a minimum miss distance of 0.003344 at a time 
ol 2.63 sec. 

2. Adjoint Solution for Second Order System 

In the adjoint domain the system will travel in negative time starting at 
the origin and traveling outward to the initial conditions for the forward time 
system. Because this is a time invariant system the trajectory for the adjoint 
solution will be the same as the forward time system but in the opposite 
direction. From our second order example, the optimum switching in negative 


time from the terminal state at the origin iS 


(аъ) 


t-t. = 5.13 
f S N ( ) 


The trajectory constraint yields 


x (ts ха (| 


2N 


х (1) == 


A reverse impulse of twice the magnitude is applied and the trajectory 


(5.14) 


proceeds out to the desired initial conditions given. The negative time impulse 


response is thus prescribed by using the switching times 


EN = — |N x (0) +4 x7 (0) (5.15) 
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мита = МУ x, (0) 5 x, (0) , (5.16) 


The parameters for the adjoint solution are the same as for the 
forward time, but the initial values for the states are at the origin. The output 
of the adjoint system, according to (5.12), is ps, so the phase plot of Figure (5.8) 
plots -p? and ps. The adjoint solution traces almost the same path as the 
forward time solution, switching at t, 2 1.582 sec., and missing the point (2,1) by 
0.002733 at tp = 2.83 seconds. By changing the initial sign of the control effort, 
and adjusting the switching time, we could drive the adjoint system to any 
desired point in the phase plane, and therefore know the optimal solution for the 


forward time system. 
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Forward Time Phase plot 


min miss = 0.006596 
ts = 2.581 sec. 
tf = 4.165 sec. 





Figure 5.7 Simulation of Forward Time System 


Adjoint Solution Phase 


min miss = 0.002733 
ts = 1.582 sec. 
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Figure 5.8 Simulation of the Adjoint Solution 


50 


D. MISSILE SIMULATION WITH ADJOINTS 

The method of adjoints can be extremely useful when dealing with time 
varying systems. By using the adjoint solution we are able to see how the 
forward time system behaves for all times. 

1. Missile/Target Model 

We will simplify our scenario of missile/target engagement in order to 

present an uncluttered example of the method. The target will have zero x- 
velocity and constant y-velocity, and will start on the x-axis at a distance, R, 
from the origin (see Figure 5.9). The missile will start at the origin with a 


constant x-velocity by zero initial y-velocity. 


Missile 





Figure 5.9 Geometry of Missile/Target Adjoint Solution 
Since R is large and y,- y, is small we use the small angle 
approximation and define 
o = )م‎ ¥ 0, (2:095 
R R 
The closing velocity 1s constant so we note that 


R = V,(t, - t) (5.18) 


ЭЛ 


where tr 1s the length of the simulation. 
2. Signal Flow Diagram 
Because the x-velocities are constant we need only to model the y- 
dimension. We develop the flow diagram in Figure 5.10 from the geometry of 
Figure 5.9. We use proportional navigation control and impulsive inputs are 


used for initial conditions. 


O(t) ГА p? pA Estimator 


Or 
Observer 





Figure 5.10 Signal Flow Diagram for Missile/Target Model 
The estimator used was developed from the mechanics and dynamics 
of a seeker head system, however, it can be shown to be equivalent to a 
Kalman Estimator or Luenberger Observer. This estimator in Figure 5.11 has a 
time constant of 0.1 seconds. The primary interest in this model is to determine 


the miss distance at the final time, t = tr, so the output of the system IS ymiss- 


J2 





Figure 5.11 Signal Flow Diagram of the Estimator 


3. Forward Time Simulation 


We define our states to be 


А] Ут 
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машна 000) 


Since the x-velocities are constant we know approximately when the 
point of closest approach will occur, and therefore about how long to run the 
simulation. The LOS angle, o, is a function of t and tr, so that different values of 
tr will result in changes to the control effort and system response. Because this 
state matrix, A, is time dependent , it must be redefined at each time step of the 
simulation. We define our proportionality constant, n, to 4, the closing velocity, 
Ve, to 5000 ft/sec., so that with tr = 4 sec. the initial range is 20,000 ft. The 
output of the simulation is shown in Figure 5.12. 

To study the effect of different values of tr, or to locate the optimal 
value, we may have to run the forward time simulation many times, which 
could be a very tedious process, especially since we are only interested in the 
final value of t = tp. 

4. Adjoint Solution 

An alternative to running the forward time system over and over again 
would be to run the adjoint solution once to generate the final values of the 
family of forward time solutions. This time we generate the adjoint system 
using matrix algebra instead of the graphical method. The solutions to the 


adjoint system are, from (5.9) and (5.10) 
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—100 








0 0 О О О 
t= 
l 0 0 0 0 0 -] 
0 
0 0 0 —100 0 0 0 
pee: ko; O iM 
= 
0 0 0 0 l 0 
ОООО и (5.23) 


This system generates a curve whose value at any time, t,, is the final 
value of the forward time system where ts =t,. Figure 5.13 shows four curves 
from the forward time solution corresponding to ts = 1,2,3, and 4 sec. Each 
curve end exactly on the adjoint solution curve at the corresponding adjoint 


time. 
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Figure 5.12 Forward Time Simulation 
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Figure 5.13 Adjoint Solution Curve 
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VI. CONCLUSIONS AND COMMENTS 


We have developed the second order minimum time optimal controller and 
applied it to the fast reaction missile defense problem with both the closed form 
analytic solution, and the open form switching curves solution. In comparing 
our open form solution with a standard proportional navigation controller 
solution we demonstrate an increased maneuverability enabling our missiles to 
defend against faster, more capable threats. The accuracy at intercept 1s a 
function of the control logic used to shut off the control effort when desired 
conditions are met. This is a subject that should be explored in future projects. 

We introduced a third order minimum time controller which promises to not 
only improve reaction time and maneuverability, but also presents us with the 
ability to control and define a desired attack angle for an improved destructive 
potential. A model for a practical third order controller should be developed 
and evaluated. 

Having introduced the method of adjoints and shown some of its functions, 
we hope to stimulate some more practical applications of this technique. We 
are excited at the possibilities of using the method of adjoints in determining 
optimal, closed form solutions to forward time problems at speeds fast enough 


for practical implementation. 
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APPENDIX A. PROGRAM CODE 


All of the simulations for this project were run on both IBM-AT! class and 
Macintosh II? computers using the matrix manipulation language MATLAB. 
For IBM-AT compatibles MATLAB, version 3.5f was used, and on the 
Macintosh II computers version MAC II-MATLADB, version 1.1b. This 
appendix contains the source code for all of the simulations and functions 
written in support of this project. 

Only a limited background in programming is required for understanding 
these files. While MATLAB is similar to FORTRAN, MATLAB's control 
Structures are much less complex, and with matrix manipulation built into the 
system, vector definition and storage are greatly simplified. Comments are 
Started by the percent sign (%) and continue to the end of the line. Ellipsis (...) 
at the end of a line indicate the continuance of the logical line onto the next line 
of code. 

Each file will begin on a new page to assist those who are interested in 
examining or reproducing the code. Although an analysis of these files is not 
necessary to understand this report, readers are encouraged to examine them 
closely for further information. 

The code is presented in the order of usage in the main text with all 


supporting functions grouped with the main program of interest. 


1 IBM and IBM-AT are registered trademarks of IBM. 
? MAC II is a registered trademark of APPLE. 
3 MATLAB isa registered trademark of The MathWorks, Inc. [4] 
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1. BB2NDST.M 


% BB2NDST.M 11 Mar. 1991 
% BB2NDST.M 1s a simulation of the 2nd order bang-bang controller with an analytic 
% solution for the Switching Time of the control effort. The signs of the control effort 


% must be matched with the initial conditions to have convergence. 


% written by Colin R. Cooper 


% Define the State Equations. 


A = [0 1:0 0]; 

В =[01]; 

ej 

x10 2- O6 x] initial condition. 
x20 = 1; % x2 initial condition. 


N =1; % Maximum control effort. 


aN = abs(N); 

Tf = 4.17; % Length of simulation. 

dt = .005; 0o Time increment for simulation. 
[Phi,Del] = c2d(A,B,dt); % Discretized System. 


% Create the storage vectors. 
kmax = Tf/dt + 1; 

x = zeros(2,kmax); 

y = zeros(1,kmax); 

time = zeros(1,kmax); 


MEM = [x10:x201; % Set initial conditions for x. 


% Define the Switching Time for the control effort. 
tsf = abs(x20)/aN + sqrt(abs(x10 + x20*abs(x20)/(2*aN))/aN); 


% Begin simulation loop. 
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[ог 1  Jikinax- ] 


if time(1)«tsf 
u N. 
else 
u=N; 
end 


x(:,i1+1) = Phi*x(:,1) + Del*u; 
УТЛЫ: 
time(i+1) = time(i)+dt; 


end 


% Plot the output of the simulation. 


clg,plot(x(1,:),x(2,:),'-w’'), grid 
title('Phase plot - Switching Time’) 
xlabel(X1'),ylabel(X2') 
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2. BB2NDSL.M 


% BB2NDSL.M 11 Mar. 1991 
% is a simulation of the 2nd order bang-bang controller using a Switching Law for the 
% control effort. 


% written by Colin R. Cooper 


A = [0 1:0 0];% Define the State Equations. 

E OU; 

STO 

xO = 2: 9o x1 initial condition. 

к=]: % X2 initial condition. 

N = 1; % Maximum control effort. 

Tf = 4.5; % Length of simulation. 

di = ОГ % Time increment for simulation. 
[Phi,Del] = c2d(A,B,dt); % Discretize the System. 


% Create storage vectors. 

kmax = Tf/dt + 1; 

x = zeros(2,kmax); 

y = zeros(1,kmax); 

u = zeros(1,kmax); 

meee) = (x 10:x 201: % Set initial conditions for x. 


% Begin simulation loop. 

for1 = 1:kmax - 1 
u(i) = -N*sign(x(1,1) + .5*x(2,1)*abs(x(2,1))/N); 
x(:,1+1) 2 Phi*x(:,1) + Del*u(i); 
(1441) 2 C*x(:,14-1); 


end 
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% Plotthe output of the simulation. 
clg,plot(x(1,:),x(2,:)),grid 

title( Phase plot - Switching Law’) 
xlabel( X 1),ylabel( X2") 
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3. 


% 


% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 


% 


SIM2SL.M 
SIM2SL.M 2nd Order Switching Laws Control 08 Apr. 199] 
Simulation of the missile/target simulation. 


-> Uses analytic values for sigma-dot,sigma-ddot for calcultion of the control effort 
using the 2nd order switching curves. 

-» The Target makes no evasive maneuvers. 

-» Calculates the Quantization Error based on the average velocities for the crossover 
endpoints. 

-> Calls INTERP.M function which takes the states at the crossover endpoints, creates 
100 point straight lines to connect the points, and determine the minimum miss 
distance. 

-> Allows delay time for missile control. 


-> Target is now on level flight with Beta = 0. 


written by Colin R. Cooper 


7o Define states. 

N = 1000; 9o Maximum control effort. 
Am = [01 0 0;0 0 0 0;0 0 0 1;0 0 0 0]; 7o Missile State Equations. 
Bm = [0 0;1 0;0 0;0 1]; 

At = [01 0 0;0 00 0;0 00 1;00 0 0]; % Target State Equations. 
Ji C | 001000 01]: 

mnm 25. % Total time of simulation. 
dps 201; % Sample step size. 

7o Descretize the states. 

[Phim,Delm] = c2d(Am,Bm,dt); 7o Discrete Missile System 
[Phit,Delt} = c2d(At,Bt,dt); % Discrete Target System 
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% Define storage vectors. 
kmax = Tf/dt + 1; 
xm = zeros(4,kmax); % xm=[y yd x ха |] 


xt = zeros(4,kmax); % xt=[y yd x xd ]' 


time = zeros(1,kmax); 


sig = zeros(3,kmax); 


am = zeros(1,kmax); 


um = zeros(2,kmax); 


Rim = zeros(1.kmiax): 


xm(:,1) = [0;0;0;2000]; % Initial Conditions for Missile. 
xt(:,1) 2 [1000;0;10000;-3000]; % Initial Conditions for Target. 
Run(1) 2 sqrt((xt(1,1)-xm(1,1))^2 + (xt(3,1)-xm(@,1))42); 9% First value for Range. 
acm = 0.0; Zo Initial acceleration for the missile. 
act = 0 0 Zo Initial acceleration for the Target. 


7o Begin the Simulation Loop. 


for1= l:kmax-1 


% Define angles. 


Jo 


beta = atan2(xt(2,1), xm(4,1)); % Velocity angle for the Target. 

gam - atan2(xm(2,1), xm(4,1)); 7o Velocity angle for the Missile. 

sig(1,i) = atan2(xt(1,i)-xm(1,i), xt(3,1)-xm@,i)); 96 LOS angle. 

S1g(2,1) 2 ((xt(3,i1)-xm(3,1))* (xt(2,1)-xm(2,1)) - (xt(1,i1)-xm(1,1))... 
*(xt(4,1)-xm(4,1)))/Rtm(i)‘2; 

S1g(3,1) 2 ((xt(3,1))-xm(3,1)) *(act*cos(beta)-acm*cos(gam))-... 
(xt(1,1)-xm(1,1))*(act*sin(beta)+acm*sin(gam))+... 
2* (xt(4,1)-xm(4,1))*(xt(2,1)-xm(2,1)))/Rtm(@)42; 

sig(3,1) = sig(3,1)+2*(xt(4,1)-xm(4,1)+xt(2,1)-xm(2,1))*sig(2,1)/Rtm(); 


Acceleration = Max value perpendicular to gamma. 
am(1) = N*sign(sig(2,1) + sig(3,i)*abs(sig(3,1))/(2*N)); 
um(:,1) = [am(i)*cos(gam); -am(i)*sin(gam)]; 
xm(5i-l1)-PhmsxmO Delum 1: 

xt(:,1- 1) = Phit*xt(:,1) * Delt*[0;0]; 
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time(i+1) = time(i) + dt; 
Rtm(i+1) = sqrt((xt(1,i+1)-xm(1,i+1))A2 + (xt(3,1i+1)-xm(3,i+1))A2); 
acm = am(1); 


end 


% Evaluation of Miss distance. 


iflag = 0; 

rm = find(Rtm==min(Rtm)); % Index of the minimum range value. 

if rm == kmax 7o Determine whether the crossover occurs 
iflag = 1; % before or after the min range value. 
It=rm- 1; 


elseif Rtm(rm+1)<Rtm(rm-1) 
It = rm; 

else 
It = rm - 1; 


end 


% Average Velocities in Intercept Area. 

Vm z .5*(sqrt(xm(2,1t)^2--xm(4,It)^2) 4 sqrt(xm(2,It- 1)^24xm(4,It- 1)^2)); 

IE" (sqrtOu2. IU^2-* x2: 1^2) 4 sqrt(xt2,It- D^23*xt(4, It D)^2)); 

QE = .5*dt*(Vm + Vt); % Quantization Error. 

if iflag == 
г = interp(xm(:,It:It+1),xt(:,It:It+1)); % Interpolated minimum miss distance. 
ISG = [Inter Miss= num2str(r) fr]; 

else 
rstr = ['Crossover never reached’; 

end 


% Plot the output of the simulation with results. 
plot(xm(3,:),xm(1,:),'-w’,xt(3,:),xt(1,:),'-w’) 
text(.15,.85,[l'Intercept Time = ' num2str(time(rm)) ' sec’],’sc’); 
text(.15,.81,['Min. Miss 2 ' num2str(Rtm(rm)) ' ft'],'sc); 
fext(-15,,/7.| Quan. Err = “num2str(QE) ft |, se); 

TEX 15. о со SCY 
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(е5 Со = пише о па so); 
title('Control: 2nd Order Switching Laws (sig, sigd)') 
xlabel(Direction 1 (ft)),vlabel( Direction 2 (ft)) 
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4. SIMPRNV.M 


% SIMPRNV.M simulation of the missile/target simulation. 


% Proportional Navigation Guidance. 


% written by Colin R. Cooper 


% Define states. 

As = [0 1; -64 -16]; 

Bs = [0; 64]; 

Am =[0100;0000;0 0 0 1;00 0 0]; 
Bm = [0 0;1 0;0 0;0 1]; 

At=(0 100;0000;000 1;000 0}; 
IE 00:1 0:00:0 LT 

t= 2.25; 

duas: 

maxac = 1000.0; 


% Descretize the states. 
[Phis,Dels] 2 c2d(As,Bs,dt); 
[Phim,Delm] = c2d(Am,Bm,dt); 
[Phit,Delt] = c2d(At,Bt,dt); 


% Create strage vectors. 
kmax = Tf/dt + 1; 

xm = zeros(4,kmax); 

xt = zeros(4,kmax); 
um=zeros(2,kmax); 
time = zeros(1,kmax); 
Кип = zeros(1,kmax); 


p = zeros(2 kmax): 


% Estimator for sigma-dot. 


% Missile State Equations. 


% Target State Equations. 


% Simulation Time. 


% Sample step size. 


% Estimator System. 
% Missile System. 
% Target System. 


% xm = [ y yd x xd | 
% x= [y vd x xd | 
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xm(:,1) = [0:0:0;2000]:; % Missile Initial Conditions. 
xt(:,1) = [1000;0; 10000;-3000]; % Target Initial Conditions. 
Rtm(1) = sqrt((xt(1,1)-xm(1,1))42 + (xt(3,1)-xm(3,1))^2); 96 First value for Range. 


% Begin simulation loop. 
fori = l:kmax-1 
Vm = sqrt(xm(2,1)42 + xm(4,1)*2); 
ub(1) = atan2(xt(1,1)-xm(1,1), xt(3,1)-xm(3,1)); 
b(:,i+1) = Phis*b(:,1) + Dels*ub(i); 
gam - atan2(xm(2,1), xm(4,1)); 
if Vm*b(2,1)<=maxac % Apply Thrust limitation. 
an = Vi bo Dp: 
else 
ат = тахас; 


end 


um(:,1) = [am*cos(gam); -am*sin(gam)); o Acceleration is perpendicular to gamma. 


xm(:,1+1) = Phim*xm(:,1) + Delm*um(:,1); 

xt(:,1+1) = Phit*xt(:,1) + Delt*[0;0]; 

time(i+1) = ttme(i) + dt; 

Rtm(i+1) = sqrt((xt(1,i+1)-xm(1,1+1))42 + (xt(3,i+1)-xm@3,i+1))%2); 


end 


rm = find(Rtm==min(Rtm)); % Index of minimum range value. 
if Rtm(rm-* 1)«Rtm(rm-1) o Determin whether crossover occurs 
[t = ru % before or after the min range value. 
else 
К =пп- |; 
end 
г = interp(xm(;, It:It+1),xt¢,It:It+1)); % Obtain the Interpolated min miss distance 


Vm = sart(xm(2,It)42 4 xm(4,It)^2); 
Vt = sqrt(xt(2,It)A2 + xt(4,It)A2); 
QE = .5*dt*sqrt( Vm^2  Vt^2); % Quantization Error. 


7o Plot the output of the simulation. 


axis([O 10000 0 1400]); 95 Same scale as the Optimal Control Sim. 
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POMO, xml) ow XIO, XII, W) 
text(.15,.85,('Intercept Time 2 ' num2str(time(rm)) ' sec'],'sc?) 
text(.15,.81,['Miss distance = ' num2str(Rtm(rm)) ' ft'],'sc'); 
text(.15,.77,l'Quant Error 2 ' num2str(QE) ' ft'],'sc); 
text(.15,.73,[ Inter Error = ' num2str(r) ' ft'],'sc'); 
text(.15,.69,['Sigma(It) 2 ' num2str(ub(It)* 180/pi) ' deg'],'sc'); 
text C 15,.65, [Max Acc. Limit = ' num2str(maxac) ‘ ft/sec“2’], sc’); 
title('Control: Proportional Navigation") 

xlabel(' Direction 1 (ft)),ylabel( Direction 2 (ft) ); 


69 


5. BB3RDORD.M 


% BB3RDORD.M 3/7/91 
7o This is a simulation of a 3rd order system, forward time. 


Ço This version allows varied N values. 


% written by Colin R. Cooper 


х10=-.5; % Define initial conditions 

x20=0; 

x30=0; 

N=1; 9% Set the magnitude of the control effort. 

Tf=2.5; Zo Set maximum time of simulation. 
dt=.002; 2o Set the simulation step size. 
A=[0 1 0;0 0 1;0 0 0]; % Define the State Equations 

В=[0 0 1]; 

CET OO 

(Phi,Del]=c2d(A,B,dt); 2o Discretize the system. 
kmax=Tf/dt +1; % Max integer value for the simulation. 


7o Prepare storage vectors. 
u=zeros(1,kmax); 
X=ZEIOS(o ex) 
y=zeros(1,kmax); 
w=zeros(1,kmax); 
f=zeros(1,kmax); 


time-zeros(1,kmax); 
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ОТО хоо хз % Define initial conditions in state vectors 


% Begin loop for simulation. 

for (121:kmax-1) 
w(i)ssign(x(2,1)*x(3,1) *abs(x(3,1))/(2*N)); % Defining the switching law. 
f(i)=x(2,1)+w(i)*(x(3,1)42)/(2*N); 
u(1,1)=-N*sign(x(1,1) + (x(3,1)43)/(3*N42) + wG)*x(3,1)*x(2,D/N +... 


f(1)*abs(f)/N)^.5); 
X(:,1+1) = Phi*x(:,1) +Del*u(1,1); % Calculate the state values. 
y(1,i+1) 23 C*x(:1*1); 
time(141)2 time(1) + dt; 9o Store time vector. 


end; 


2o Plot the switching law and its components. 

clg, axis([O 2.5 -1.5 1.5]); 

plot(time,u,time,.75*w,time,f) Ç w is scaled to distinquish it from u. 
title(Control Laws vs Time") 


pause 


2o Plot the 3-Dimensional view of the simulation from 45 deg. azimuth 
o and 45 deg. elevation angle. (Pos. x1 vector is out of the screen 

% and towards the left, Pos. x2 is out of the screen and towards right. 
clg, axis; 

РО ОСА) ХО) ХК) 45,45): 
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6. ADJOINT.M 


% ADJOINT.M is a simulation of a controlled system using formard time simulation 


% and reverse time or adjoint solution. 


% written by Colin R. Cooper 


A = [0 1 0;0 0 1;0 0 0J; 


В =[00 1]; 

С=[100]; 

х10 = 2: 

x20 2 0; 

М = 1; ?c Maximum control effort. 
І = 4; 

аї = .005; 


tsf = x20/N + sqrt(x10/N + .5*(x20/N)^2); 
tsa = sqrt(x10/N + .5*(x20/N)^2); 
kmax = Tf/dt + 1; 


% Define the Storage Vectors 
X = zeros(3,kmax); 

y = zeros(1,kmax); 

time = zeros(l,kmax); 

imp - zeros(1,kmax); 

imp(1) = -N/dt; 
imp(round(tsf/dt)+1) = 2*N/dt; 
x COND [x 10-20); 
[Phi,del]=c2d(A,B,dt); 


for 1 = 1:kmax - 1 
x(:,1+1) = Phi*x(:,1) + del*imp(i); 
y(1,+1) = C*xC,i*1) 


22 


11 Mar. 199] 


7o Define Forward Time State Equations 


% x1 initial condition. 
00 х2 initial condition. 


% Length of simulation. 

7o Time increment for simulation. 
% Forward Time Switching Time 
% Adjoint System Switching Time 
% Max length of storage vectors. 


do States of the system 
% Output state vector 


% Impulse vector for control times. 
% First pulse at t = 0 

% Second pulse at t = tsf 

Ç Set initial conditions for x. 


2o Discretize the state equations 


% Begin the simulation loop 


time(i*1) 2 time(1)«dt; 
end 
miss = mnsan ха |) 2) 


tff = time(find(sqrt(x(1,:).42 + x(2,:).42)== ... 


Imin(sart(x(1,:).A2 + x(2;:).A2)))): 
clg,plot(x(1,:),x(2,:),'-w'), grid 
title((Forward Time Phase plot) 
xlabel( X1 ),ylabel( X2") 


text(.6,.80,['min miss 2 ' num2str(miss)],'sc') 
ruo ЕКЕ УО ВЕС SC) 


C674 (UE = ооа ГЕ SEC SC) 


pause 


7o Now for the Adjoint System. 
xa = zeros(3,kmax); 

time - zeros(1,kmax); 

impa = zeros(1,kmax); 

impa(1) = N/dt; 
impa(round(tsa/dt)+1) = -2*N/dt; 
[Phi,del]2c2d(A',C',dt); 


јог1 = 1:kmax - | 
xa(:,1+1) = Phi*xa(:,1) + del*impa(i); 
ya(1,i+1) = В'*ха(:+1); 
time(it+1) = time(a) + dt; 


end; 


missa = min(sqrt((xa(3,:)-x10).42 + (xa(2,:)-x20).42)); 


plot(xa(3,:),-xa(2,:),'-w'), grid 
title Adjoint Solution Phase plot’) 
xlabel('X3'), ylabel('X2') 


Йо 


% Find the min. miss distance 


96 Find the time of the min. miss 
9o Plot the Phase Plane 


% Lable the graph and display the desired 


% Information 


% Define storage vectors 


% First impulse 
% Second impulse 


% Discretize the Adjoint System 


% Begin the simulation loop 


% from the initial conditions 
tfa = time(find(sart((xa(3,:)-x10).42 + (xa(2,:)-x20).42)== ... 
min(sqrt((xa(3,:)-x10).^2 + (ха(2,:)-х20).^2)))); % 


9o Find the time of the 


min miss distance 


7o Plot the Phase Plane 
2o Label the graph and display the desired 


% Information. 


% Find the min. miss distance 


text(.6,.8,['min miss = ' num2str(missa)],' sc") 
text .6..77 tsa = оа) sec se) 
text(.6,. 74, [tia mum2sa (Ga sec Se 
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7. ADJSIM.M 


% ADJSIM.M Adjoint solution to missile intercept problem. 


% written by Colin R. Cooper 


n4 

У = 5000; 

к= 1000001]; 
cC —-|[- 100010] 
Her 

Gi Ol 


11 Mar. 1991 


2o Proportionality Constant 
% Closing Velocity 
% Define State Matrices 


% Set first value for incremented times 


% This outer loop runs the complete forward time system for each value of Tf, storing 


p the output into vectors at the end of the loop. 


for j = 1:4 
kmax = Tf/dt + 1; 
imp = zeros(1,kmax); 
imp(1) = 1/dt; 
= zeros(6,kmax); 
y — zeros(l,kmax); 


ПС = zeros( 1 kmax): 


% Forward Time Simulation. 

count = 0; 

(ог1 = 1Кктах - 1 
count = count + 1; 
ki = 1/(v*(Tf - итед) + 1е-12)); 
A=[010000 
00n*v 000 
000100 


D 


Zo Maximum index for storage vectors 
% Create vector for impulsive input 
% Define the pulse at t = Û 


% Create storage vectors for states 


% Counter to indicate the computer is busy 


% Begin simulation loop 


% 1/Range 
% Define the time varying A matrix 


-100*ki 0 -100 -20 100*ki 0 
000001 
000000); 


Phi = eye(6) + A*¥dt + .S*AA2*dtA2; 


x(:,1+1) = Phi*x(:,) + B*dt*imp0(i); 
y(:,1+1) = C*x(:,1+1); 

time(i+1) = time(i) + dt; 

lÍ count == 90 


disp( working; 
count = 0: 
end 
end 
Тела) = ТЕ; 


yend()) = y(length(y)); 
eval y nuim2sir(jp = у); 
eval([t',num2str(J), = time;']);1 
Tf = Tf+ 1 


end 


% Adjoint Simulation. 

J = 4 5 

kmax = Tf/dt + 1; 

imp = zeros(1,kmax); 

imp(1) = 1/dt; 

ха = zeros(6,kmax); 

уа = zeros(1,kmax); 

time = zeros(1,kmax); 

count = 0; 

fori = l kmax i 
count = count + 1; 
ki = 1/(v*(time(i) + 1e-12)); 
A=|010000 
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% Discretize the matrix with 2nd order 


o expansion - drop higher terms 


2o Counter indicates computer is busy 


% End simulation loop 


2o Record final time 

% Record corresponding final value ymiss 
% Save output vector 

% Save corresponding time vector 

% Increment to next Tf 


2o End outer forward time loop 


% Set a simulation time 

o Maximum index for vectors 

Ço Create vector for impulsive input 
% Define the impulse at t = Û 


Ç Create storage vectors 


% Begin simulation loop 


% 1/Range 
% Define the time varying A matrix 


00n*v000 
000100 
-100*ki 0 -100 -20 100*ki 0 
000001 
000000] 
Phi 2 eye(6) * А" “а + .S*A'^2*gt^2; 7o Discretize the matrix 
ха(:1+1) = Phi*xa(:1) + C'*dt*imp(); 
ya(:,1+1) = B'*xa(:,1+1); 
time(i+1) = time(i) + dt; 
if count == 50 
dispC working’); 
count = Û; 
end 
end % End adjoint simulation loop 


plot(tl,yl,'-w',t2,y2,'-w',t3,y3,'-w ,t4,y4,'-w’,time,ya,'-w’) % plot output 
xlabel(‘Time (sec)'),ylabel( Ymiss (ft)'), grid 


E 


8. INTERP.M 


function r = interp(xm,xt) 
% INTERP.M will return an interpolated value for the min.miss distance given the states 


% for the two intercept values. 


7o Written by Colin R. Cooper 26 Mar 1991 


% Increase the sample rate by a factor of 100 in crossover region. 
dax 2 (xm(3,2)-xm(3,1))/100; 

day = (xm(1,2)-xm(1,1))/100; 

dbx = (xt(3,2)-xt(3,1))/100; 

dby = (xt(1,2)-xt(1,1))/100; 


% Define the storage vectors for 100 data points. 
ax = zeros(1,100); 
ay = zeros(1,100); 
bx = zeros(1,100); 
by = zeros(1,100); 


% Set initial values for each vector. 
ax l) = XML), 

ay(1) = xm(1,1); 

рО xt T). 

by(1) 2 xt(1,1); 


% Assuming a straight line trajectory with constant velocity in the crossover region, create 
% the interpolation data sets. 
for 1 = 1:99 

ах(1+1) = ах(1) + аах; 

ау(1+1) = ay(i) + day; 

bx(-*1) 2 bx(1) * dbx; 
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ру(1+1) = by(1) * dby; 
end 


% Find the closest point of approach of the two line segments. 
r - min(sqrt((ax - bx).^2 + (ay - by).42)); 


uo 


9. PLOT3D.M 


% PLOT3D isa 3d plotting function allowing rotation and elevation adjustments. The plot 
% shows a 3-D curve and its projection onto the X-Y Plane. 

9% X,Y,and Z data must be passed, and if no azimuth or elevation values are passed 
% they default to AZ = 45°, EL = 30°. The Azimuth is the angle of rotation of the view 

7o angle about the Z-Axis. AZ = 0° is looking straight down the X-Axis at the Y-Z Plane. 
% The elevation is the angle from which the plot is viewed. EL = 90° is looking down the 
% Z-Axis at the X-Y Plane. The elevation angle can vary from -90° to 90°. 

% The tick marks on the axis will default to 10 marks per axis and the values will be 
% displayed on the screen. To define the values of the tick marks a three element vector 
o must be passed containing [dx dy dz]; 

% Axis values will be calculated and fixed unless an axis vector 1s passed to the 

o program: [xmin xmax ymin ymax]. 

% The Transformed 2-D vectors V and Vs are returned, where V is the 3-D curve and 
o Msisthe Projection on the X-Y Plane. 

% 

% Example : [V,Vs] = plot3d(x,y,z,-45,30,dx,Ax) 


% written by Colin Cooper 4/26/91 
function [V,Vs]=plot3d(x,y,z,az,el,dx,Ax) 


if nargin < 5, el = 30; end 
if nargin « 4, az 2 45; end 


az=-az*pi/180; el=el*pi/180; 
alpha = 45*pi/180; beta=30*pi/180; % Angles for mapping onto 2D 
% alpha = rot. Beta = elv. 
S1=[sin(az) cos(az) 0 
-sin(el)*cos(az) sin(el)*sin(az) cos(el)]; 
S2-[sin(az) cos(az) 0 
-sin(el)*cos(az) sin(el)*sin(az) 0]; % Trans for Projection. 
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[r C] size(x); 
MT i 
KSA y y A. 


end 


У = 51*%[х;у;2]; 
Nu 21у 


if nargin « 6 
dx(1) = (max(x)-min(x))/10; 
dx(2) = (max(y)-min(y))/10; 
dx(3) = (max(z)-min(z))/10; 


end 


ax 1 = [min(x) max(x); 0 0; 0 0]; 
ах2 = [0 0; min(y) max(y); 0 0]; 
ахз = [0 0; 0 0; пиуп(2) тах (2) |; 


ax1t 2 [fliplr(0:-dx(1):min(x)) 0:dx(1):max(x) 
zeros([0:-dx(1):min(x) O:dx(1):max(x)]) 
zeros([0:-dx(1):min(x) O:dx(1):max(x)]p]; 

ax2t = [zeros([0:-dx(2):min(y) 0:dx(2):max(y)]) 
fliplr(0:-dx(2):min(y)) 0:dx(2):max(y) 
zeros([0:-dx(2):min(y) 0:dx (2):max(y)])]: 

ax3t 2 [zeros([0:-dx 3):min(z) O0:dx (3):max(z)]) 
zeros([0:-dx(3):min(z) 0:dx(3):max(z)]) 
fliplr(0:-dx(3):min(z)) 0:dx(3):max(2)]; 


axe 51*ах[; 
ах = 5] "ах2 
ас = 17 ax3: 
ax lit У ах и: 
ЕО ИС 
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ах s tax f: 


прах = пио (У) УУ ах ахо) ахо j; 
if minx == 0, minx = -1; end 

maxx = max([ V(1;:) Vs(1.:) ax1(1 2 ax2(15 ax (15) 
if maxx == 0, maxx = k end 

miny = min([V(2,:) Vs(2,:) ax1(2,:) ax2(2,:) ax3(2,:)]); 
if miny == 0, miny = -1; end 

maxy = max([V(2,:) Ys(2,:) ax1(2,:) ax2(2,:) ax3(2,:)]); 


if maxy == 0, maxy = 1; end 


clg 
axis('square ); 


if nargin « 7, 


axis([minx-.3*abs(minx) maxx4.3*abs(maxx) miny-.3*abs(miny) ... 


maxy+.3*abs(maxy)]); 
else 
axis(Ax); 
end 
b = axis; 
hold on 
роках ахо ах) ахо s 
ах2(1.:),ахо (о) WV QASL. VOLO TT 
ax3(1,:5ax902, 3), -w ak t(1 ахо у 
lv = length(V); 
plot(V (1,:).V @,:), -w Vs(1;:).VsO wM 
(МО ОЛИ УВ БИМ РУС АПИЈА Vs Q TENE) 
plot([b(1) b(2) b(2) b(1) b(1)],[b(3) b(3) b(4) b(4) b(Š)], -w ); 
hold off 
text(.22,.08, [AZ = ' num2strCaz* 180/pi) ° | sc > 
text(.72 08 [EL = Ио аеро Ее 
text(.22 88. [dx = пий d(I se» 
text(.30 38 Fay Жашта) se). 
кех 72. 68 рас = ооа о se 
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